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ABSTRACT 



We report observational upper limits on the mass-energy of the cosmological 
gravitational- wave background, from limits on proper motions of quasars. 
Gravitational waves with periods longer than the time span of observations 
produce a simple pattern of apparent proper motions over the sky, composed 
primarily of second-order transverse vector spherical harmonics. A fit of such 
harmonics to measured motions yields a 95%-confidence limit on the mass-energy 
of gravitational waves with frequencies v < 2 x 1CT 9 Hz, of < 0.11/i~ 2 times the 
closure density of the universe. 



Subject headings: Cosmology: Gravitational Radiation, Gravitation - 
Techniques: Interferometry 
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1. INTRODUCTION 

Like variations in refractive index, the changes in spacetime metric produced by 
gravitational waves alter optical path lengths. Light propagation through gravitational 
waves preserves sources' surface brightness and total intensity, to first order in the wave 
amplitude ( Pondi, Pirani, fc Robinson 1959| , [Zipoy 1966| , Penrose 1966| ); but it can produce 



oscillations in apparent position, at the period of the gravitational wave QFakir 1994| , purrer 



1994| , |Bar-Kana 1994| , |Pyne et al. 1996| , |Kaiser fc Jafle 1996j ). Over time intervals much 



shorter then a gravitational wave period, these deflections cause a characteristic pattern of 
apparent motions in the plane of the sky (in astronomical parlance, proper motions) ( |Pyne 



|et al. 1996|) . In this Letter we set stringent upper limits on the energy density of such 
waves, using measurements of the proper motions of quasars. 

Although detected only indirectly to date, most cosmologists believe gravitational waves 
are commonplace. Very low-frequency gravitational waves arise naturally in inflationary 
cosmologies ( [Rubakov fc Sazhin, fc Veryashkin 1982| , |Fabbri fc Pollock 1983j ) and dominate 



the mass-energy of the universe in some of them QGrishchuk 1993|) . Other possible sources 



of such a background include phase transitions in the early universe, networks of cosmic 
strings, and collisions of bubbles QVilcnkin 1981 , Hogan 1986 ). 



Various observations set direct or indirect observational constraints on the spectrum 
of low-frequency gravitational radiation. These constraints are often expressed in terms 
of QgW, the ratio of the mass-energy density of the gravitational waves to that required 
to close the universe. Some observational constraints are sensitive to a narrow range of 
frequencies, and are best expressed in terms of the logarithmic derivative dQcw / 'dlnu 
evaluated at some frequency v. Timing of pulsars sets observational limits for periods 
as great as the span of observations (|Backer fc Hellings 1986| ). For the extremely 



stable millisecond pulsars the limit corresponds to d£L GW jd\nv < 10 , at frequency 
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v = 4.4 x 10~ 9 Hz ( |Kaspi, Taylor, fc Ryba 1994]) . After taking into account the energy 
carried away by gravitational waves and the effects of galacto centric acceleration, orbital 
periods of binary pulsars are sensitive to gravitational waves with periods as great as 
the light travel time from the pulser. Current limits from such data sets a limit of 
tt GW h 2 < 0.04 for 10- 11 <i/<4.4x 10~ 9 Hz, and of VL GW h 2 < 0.5 for 10~ 12 < u < lO" 11 Hz 
( |Bertotti, Can, k Rees 19S5[ [l'aylor fc Weisberg Egg , fl'horsett & Dewey 199BQ . Here, 



the normalized Hubble constant is h = i?o/100 km s _1 . These measurements are sensitive 
to gravitational waves with periods as great as the light travel time from the pulsar. A 
cosmological background of gravitational waves at the epoch of recombination can produce 
anisotropies of the cosmic background radiation ( |Linder 1988a| , |Krauss fc White 1992|) . For 



gravitational- wave spectra typical of inflation, the anisotropy detected by COBE yields 
a limit of dVtow I 'd\n.v < 10 -11 , with sensitivity of the measurement concentrated near 
A ~ 2 Gpc, or v ~ 10~ 17 Hz ( Par-Kana 1994|) . Few constraints have been set in the frequecy 



range 10~ 17 < v < 10~ 12 Hz Linder (1988b) found that a gravitational-wave background 
in this frequency range would affect the galaxy-galaxy correlation function, but robust 
calculations would require knowledge of the "true" correlation function, observable only in 
the absence of such a background. 

In contrast to these techniques, proper motions of extragalactic radio sources are 
sensitive to waves of arbitrarily long wavelength, and are independent of the spectrum 
and source of gravitational-wave radiation. Under the assumption that the spectrum of 
gravitational waves is stochastic, the squared proper motion of distant sources, suitably 
averaged over the sky, is proportional to the energy density of the waves, at frequencies 
from the inverse of the period of observation to the Hubble time. 



In a previous paper (|Fyne et al. 1996] ) we discuss the pattern of apparent proper 



motions of distant radio sources, produced by a gravitational wave. We use the assumption 
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that the wavelength is short compared to the distance to the source, and the adiabatic 
approximation that the wavelength is less than the Hubble length; however, because the 
observed effect is local to the Earth, these assumptions can probably be relaxed without 
much change in the result. The best present measurements of proper motions attain 
accuracies of a few microarcseconds (//as) per year ( Eubanks et al. 1996| ) corresponding to 
observational limits on the energy density of gravitational waves of about that required to 
close the universe. 

2. THEORETICAL BACKGROUND 

Very-long baseline interferometry (VLBI) measures positions of radio sources by 
measuring the difference in arrival times of their signals at antennas in different geographical 
locations flShapiro 1976| ). The interferometrist assumes that the observations are made 
in a locally Minkowski reference frame (allowing for the orbital acceleration of the Earth 
and the general-relativistic light-bending of the Sun and planets), and so interprets these 
observations in the Gaussian normal reference frame. The delay T between arrival times 
measures the projection of the unit vector pointing toward the source onto the spacelike 
baseline vector that connects the antennas. Measurement of the delay for many sources on 
several baselines allows solution for both source positions and lengths and orientations of 
the baselines. 

Pyne et al. (1996) describe the effect of a gravitational wave on a VLB interferometer: 
the wave produces variations in delay T, which are interpreted as variations in source 
position. A gravitational wave traveling toward +z, with the "+" polarization, produces 
metric perturbations /icos(pt)(xx — yy) in the background coordinate reference frame, 
where h is the dimensionless strain of the wave, p is its angular frequency, and t is time. 
In the interferometrist 's Gaussian normal frame, the observed proper motion jl of a radio 
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source at position (9, 0) will be: 

jl = — sm(pr]) sin 9 (cos 2(p§ — sin 200) (1) 

Here, 9 measures angle from +z, the direction of propagation of the wave, and <fi measures 
azimuthal angle around it, from the x-axis; the associated unit vectors on the sky are 9 
and (f>. Proper time in the Gaussian normal frame is rj. We take h to be real, and allow 
the origin of time rj to express the phase of the wave. Fig. 1 shows the pattern of proper 
motions that this gravitational wave produces. 

The properties of this pattern of proper motions are not simple under rotation 
or superposition. However, the transverse vector spherical harmonics Y^ M ^ form an 
orthonormal basis for vector fields on a sphere, with well-understood behavior under 
rotation and superposition. Expanded in such harmonics, Eq. (|IJ) takes the form: 



jl — ph sm(pr)) 

/ 07T 

6 



x{ +^[Y 2 E +2 + Y 2 %-Y 2 M +2 + Y 2 M _ 2 ] ( 2 ) 



'[^3,+2 + Y 3j „ 2 Y 3j+2 + Y 3 _2\ + ...}. 



60 

Here we use convention of Mathews (1962,1981) for transverse vector spherical harmonics. 
These fall into 2 categories, with one family, commonly denoted "poloidal", "potential", 
or "electric" pointing down the gradients of scalar spherical harmonics; and the other, 
known as "toroidal" , "stream" , or "magnetic" pointing perpendicular to their gradients. 
We denote these categories as "E" and "M", respectively flMathews 1962| , [Mathews 198 1| ), 



with the notation (E, M) meaning "E or M" . Note that our Y^ m is the X£. m of Jackson 
(1975). The normalization condition is: 

4vr(Y £ E m Yg) = 4vr(Y £ M m Y^) = 5 u 5 mj (3) 
(YL.YJJ*) = 0. 
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Here (...) denotes an average over the sky. 

Because the Y^ M ^ form an orthonormal basis for vector fields on a sphere, they can 
represent proper motions due to a single wave, or any spectrum of waves. We use the 
expansion 



to define the coefficients m and af m in terms of the proper motion /i, observed over a 
small fraction of a gravitational wave period. The must satisfy ( [Mathews 1981|) : 



because /x is real. 



The squared proper motion averaged over the sky is related to the energy density of 
the wave Tqw by 

(a* 2 ) = y/h 2 = y (6) 

Here (..) includes an average over many wave periods, as well as over the sky. Expansion of 
the mean square proper motion in vector spherical harmonics shows that 5/6 of the mean 
square motion resides in the I = 2 harmonics: 

^Y.{KS + \<S) = % 2 )- (?) 

Because the summed, squared moduli of the a^f~ remain constant under rotation for 
for each £, and for E and M harmonics separately, Eq. (|7|) holds for a gravitational wave 
propagating in any direction with any polarization. 

If the phases, amplitudes and directions of the many superposed waves in the spectrum 
are random at the Earth, as is to be expected for a cosmological background of gravitational 
waves, the random-phase approximation shows that Eqs. (Q) and (0) will hold statistically, 
if p 2 h 2 is replaced by the integral / p 2 h(p) 2 d 3 p, and if the average of the squared proper 
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motion over time and sky is replaced by an average over an ensemble of gravitational- wave 
spectra (e.g. |Goodman 1985|) . Here h(p) 2 is the spectral density of the square of the 



dimensionless strain. Thus, Eq. (0) shows that 5/6 of the mean square proper motion due 
to the low-frequency gravitational wave background appears in second-order transverse 
spherical harmonics, and E and M harmonics contribute equally. We can express the energy 
density in such a stochastic background of gravitational waves as Slaw-, its ratio to the 
closure density of the universe, Tql = 3c 2 Hq/8tiG; and we can relate VLqw to proper motion 
in second-order transverse harmonics: 

MOW - 772^ ) - TTZ-FH 2^ 



\a 



2,m 



We seek to measure the a^^f 1 and to thus measure, or set limits on, VLqw- 



3. OBSERVATIONS AND DATA REDUCTION 

We searched for the pattern expected for gravitational waves in proper motions 
of extragalactic radio sources. We estimated second-order spherical harmonics from a 
proper-motion solution using the US Naval Observatory comprehensive database of VLBI 
observations. This contains most of the suitable astrometric and geodetic observations of 
extragalactic radio sources made with the Mark III and compatible VLBI systems |Clark et 



al. 1985|. Eubanks and Matsakis [Eubanks et al. 1996| describe the data analysis in detail 



here we present a brief overview. 

The database comprises 1,469,793 observations of delay T and its time derivative T 
with different baselines and sources. The observations were made between 3 August 1979 
and 13 February 1996, with more than 70% after 1990. Proper motions were determined for 
the 499 sources observed at more than 1 epoch. In fitting for the coefficients of transverse 
spherical harmonics we used the 323 sources known to be extragalactic by their measured 
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redshifts. Their mean redshift is 1.22. 

Reduction of this data to source positions and proper motions follows standard 
procedures for astrometric and geodetic VLBI data (e.g. Clark et al. 1985). We removed 



the best models for earth orientation, tidal deformation, and location within the solar 
system; atmospheric propagation; and solar-system effects of general and special relativity; 
before fitting for proper motion. Simultaneous observations at 2 frequencies allowed 
calibration of the propagation delay through the ionosphere. The solutions for proper 
motions also included ~ 400, 000 "nuisance" parameters, including those describing positions 
of sources and stations; clock behavior at different stations; optical path length though 
the atmosphere in a vertical direction and its horizontal gradients; unmodeled motions of 
the Earth, including variations in rotation, polar motion, nutation, and precession; and 
gravitational deflection of light by the Sun Eubanks et al. 1996| . These nuisance parameters 
are chosen to have little covariance with source and station positions, and thus to proper 
motions. After the fit, the weighted root-mean-square of the scatters of the residuals in 
delay and in delay rate were 31 ps and 93 fs s _1 , respectively. 

Effects of source structure and variations in atmospheric path length are among the 
expected systematic errors. The sources are highly energetic neighborhoods of active 
galactic nuclei. The well-known motions of jets from such sources can mimic proper 
motion, particularly in cases where the stationary core component is weak or absent (e.g. 
Guirado et al. 1995[ ). Such motions are expected to be many /ias yr _1 , but should not be 



correlated over the sky, so that they increase noise in the measurements but do not bias 
lower-order transverse spherical harmonics. Atmospheric path length is removed via models 
incorporating contemporaneous meteorological data, and the residual path length and 
its horizontal gradients are parametrized and estimated directly from the VLBI solution. 
Remaining atmospheric delay errors might show angular dependences that could mimic 
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those expected for gravitational waves. 

We then fit transverse spherical harmonics to the measured proper motions. The fit 
minimized \ 2 , the summed, squared differences of proper motions and the motions modeled 
by Eqs. (^) and (|5|). These differences were weighted by the expected errors, given by the 
standard errors from the fit to delay and delay rate, scaled by a factor of 1.35, and with an 
assumed error added in quadrature (in each component of proper motion) of 30 /ias yr _1 . 
This reweighting serves to make the \ 2 P er degree of freedom about 1, as expected; it 
accounts, in particular, for effects of variations in source structure. Comparison of fits to 
subsets of the data suggested this scaling and addition in quadrature. Fits to subsets of the 
data did not change the fitted values significantly, although the standard errors increased. 

The fit included coefficients for I = 1 transverse spherical harmonics, as well as the 
I = 2 harmonics characteristic of gravitational waves. The M I = 1 harmonics correspond 
to a rotation, of no physical significance, as it is not separable from the Earth's rotation; 
the E i = 1 harmonics correspond to acceleration not included in models used for data 
reduction. This acceleration is an interesting cosmological parameter: it is sensitive to the 
galactocentric acceleration of the Solar System, as well as any acceleration of the Milky 
Way relative to distant quasars. 

The measured values of the coefficients of second-order transverse spherical harmonics 
are statistically indistinguishable from noise. The residual x 2 to the fit was 717, distributed 
among 627 degrees of freedom. Table 1 summarizes results of the fit, and Fig. 1 shows the 
results for the i = 2 harmonics in graphical form. The fitted coefficients are combinations 
of c^m^ that satisfy Eq. (|5|) to produce real motions. The mild migration of sources away 
from the ecliptic in Fig. 1 may reflect effects of the Sun on atmospheric propagation. 



4. RESULTS 



A particular value of Slew corresponds to an ensemble of possible gravitational- wave 
spectra and coefficients a^®. Eq. (||) relates VL GW to the mean over the ensemble of the 
summed, squared coefficients. The measurement errors given in the last column of Table 1 
also contribute to the measured a^m^ • We take these contributions into account to find the 
value of VLqw that yields summed, squared coefficients less than those we observe for only 
5% of gravitational spectra, and adopt this as our upper limit. This procedure includes 
effects of sample or cosmic variance. With 95% confidence, we find that Qqw < 0.11 h~ 2 , 
where Hq = hlQQ km s" 1 . This limit holds for a stochastic spectrum of gravitational 
waves integrated over all frequencies less than half the inverse of the 10-yr span of the 
observations, or about 2 x 10~ 9 Hz. 

Considerable improvement in this limit should be possible over the next decade. The 
measurement accuracy for proper motions is proportional to the duration of observations 
to the 3/2 power. Because most data were acquired in the last 5 to 10 years, an additional 
decade of observing should improve the bound on flew by a factor of 3 to 8. Moreover, we 
can reduce effects of source structure by choosing additional sources with little structure, 
from among the thousands observable astrometrically with the Very Long Baseline Array; 
and by using models derived from images of the sources to correct for effects of structure 
directly. Together these improvements may be expected to reduce the formal errors of 
measurements of Qqw by a factor of 10 to 100, within the next 10 years. 

We thank J. Hartle for useful suggestions. This work was supported in part by the 
National Science Foundation (AST90-05038 and AST92- 17784). 
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Fig. 1.— Upper: Proper motions expected for a single gravitational wave. The metric 
perturbation is hcos(p(cz — t))(xx — yy), with z toward declination 90° and x toward 
right ascension h . Lower: Fitted coefficients of the second-order [I = 2) transverse 
spherical harmonics, displayed as proper motions at locations of sources with measured 
proper motions. Arrow lengths in degrees equal proper motion in //as yr^ 1 . Coefficients 
are not shown for £ ^ 2. Curves show the ecliptic (long dashes) and galactic plane (short 
dashes). The fitted coefficients are not statistically significant, so that the observed pattern 
of motions is consistent with filtered noise. 
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Table 1. Observed Coefficients of Transverse Harmonics 



Fitted Standard 
Parameter Value Error 

(//as vr" 1 ) (//as yr" 1 ) 



Acceleration* 

x 1.9 6.1 

y 5.4 6.2 

1 7.5 5.6 

Second-Order Transverse Vector Spherical Harmonics' 3 

a E 2 ,o 12.1 16.6 

^(a E 2 , +1 - a E 2 ,_i) -13.6 14.9 

^(aVi + «Vi) -21.7 15.6 

^(a E 2i+2 + a E 2 ,_ 2 ) 4.2 12.5 

^(a E 2 , +2 - a E 2 ,- 2 ) -0.6 13.3 

ia M 2 ,o 1.7 14.4 

-^(a M 2!+1 -a M 2,-i) 17.8 15.0 

^(a M 2 ,+i + a M 2 ,_ 1 ) -28.2 15.3 

^(a M 2 , +2 + a M 2 ,_ 2 ) -17.6 13.9 

Ma M 2,+2 ~ a M 2 ,- 2 ) 10.7 15.0 



a Effects of acceleration of the solar system barycenter relative to the observed 
extragalactic radio sources corresponds to first-order (£ = 1) harmonics. Rotations 



also correspond to first-order harmonics; our fits included these, but we do not report 
them as they cannot be separated from the Earth's rotation. 

b Second-order (£ = 2) harmonics correspond to effects of low-frequency 
gravitational waves. Otherwise identical coefficients af"^ with opposite sign of 
m have been combined to reflect the fact that the fitted motions are purely real. 



